12 Homotopy
Our final root-finding method, homotopy continuation, is a significant departure from those previously discussed. The method will ultimately require numerical integration, which will bring us into the final unit.
12.1 Basic Concepts
Homotopy describes, roughly, the process deforming one shape or object into another. In this method, we start from a complicated function and replace it with a simpler one that we can solve, then try to deform the simple problem back into the original, and follow how the root travels.
Definition 12.1 Let \(f:X \to Y\) and \(g: X \to Y\) be continuous functions. A homotopy from \(g\) to \(f\) is a continuous function \(H:X\times [a,b] \to Y\) such that \[\begin{align*} H(x,a) &= g(x)\\ H(x,b) &= f(x) \end{align*}\]
Typically, we take \(a=0,b=1\). Occasionally, we consider infinite homotopies, in which case the evaluation at \(a\) or \(b\) is replaced by a limit (in which case stronger conditions like uniform convergence could be imposed, if a pointwise limit is too weak).
We usually study functions from \(\mathbb R\to\mathbb R\) or \(\mathbb C\to\mathbb C\), in which case the following weighted averages are natural choices:
\[\begin{align*} H(x,t) &= (1-t)g(x) + tf(x)\\ H(x,t) &= (1-t^2)g(x) + t^2f(x) \end{align*}\]
Given a “difficult” function \(f\), the choice of “simple” function \(g\) is important. We can form a standard homotopy for any two functions, but we will obtain better results if \(g\) is simplified from \(f\) in a useful way. Also, we want it to be easy to find a root of \(g\). A straightforward option is \[g(x) = f(x) - f(a)\] for some fixed \(a\), since \(g(a) = f(a) - f(a) = 0\) guarantees us a known root of \(g\). Visually, this homotopy looks like sliding the graph of \(g\) up or down.
Example 12.1 Consider the “complicated” function \(g(x) = x^2 - 3\) whose root we would like to know. A “simple” but related function is \(f(x) = g(x) - g(1) = x^2 - 1\), which we know has a root \(x=1\). The standard homotopy connecting them is \[\begin{align*} H(x,t) &= t(x^2-3) + (1-t)(x^2 - 1)\\ &= x^2 - (2t+1) \end{align*}\]
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
import matplotlib.pyplot as plt
import numpy as np
from shiny.express import ui, input, render
with ui.tags.head():
ui.tags.link(
rel="stylesheet",
href="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.css"
),
ui.tags.script(src="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.js"),
ui.tags.script(src="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/contrib/auto-render.min.js"),
ui.tags.script("""
document.addEventListener('DOMContentLoaded', function() {
renderMathInElement(document.body, {
delimiters: [
{left: "$$", right: "$$", display: true},
{left: "\\[", right: "\\]", display: true},
{left: "$", right: "$", display: false},
{left: "\\(", right: "\\)", display: false}
]
});
});
""")
ui.tags.style(
"""
body { background: #C8E0FF !important; }
.irs-bar { background: #0054A9 !important; }
.irs-handle { background: #0054A9 !important; }
.irs-single { background: #0054A9 !important; color: #FFFFFF !important; }
.irs-min { background:#99CBFF !important; color: #000000 !important; }
.irs-max { background:#99CBFF !important; color: #000000 !important; }
.control-label { display: block; text-align: center; }
.shiny-input-container { margin:auto; }
"""
)
ui.input_slider(id="t",label="$t$",min=0,max=1,value=0,step=0.05)
@render.plot(alt="Graph of the slice of the homotopy at a given value of $t$.")
def homotopy_plot():
t = input.t()
def H(x,t):
return x**2 - (2*t+1)
fig, ax = plt.subplots()
xs = list(np.linspace(0.75,2,100))
ax.plot(xs,[H(x,t) for x in xs],color="#0054A9")
# formatting
fig.set_facecolor('#C8E0FF')
ax.set_facecolor('#C8E0FF')
ax.set_ylim(-1, 2)
ax.set_xlim(0.75, 2)
ax.spines.right.set_color('none')
ax.spines.top.set_color('none')
ax.spines.bottom.set_position('zero')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.legend(facecolor='#C8E0FF',frameon=False)
return figExample 12.2 You can also visualize a homotopy between univariate functions by plotting the surface it defines. Below is the standard homotopy between \[\begin{align*} f(x) &= x^3 - 3x^2 - 4x + 12\\ g(x) &= x^3 - 2x^2 - 3x \end{align*}\] with an example “slice” at \(t=0.5\).
Notice how the solutions to \(H(x,t) = 0\) form paths connecting roots of \(f\) and \(g\). The aim of the homotopy method is to find and follow this path.
12.2 Differential Equation
In our example, the solutions to \(H(x,t)\) form paths linking roots of \(f\) and \(g\). We want a way to describe these paths without directly solving \(H(x,t)=0\). To do so, we take advantage of the deformation perspective: a small increase \(\delta t\) should cause a tiny change in the graph from \(H(x,t)\) to \(H(x,t + \Delta t)\), and quantifying that kind of change is precisely what calculus does!
More precisely, one of those paths can be viewed as a function \(x(t)\) such that \(H(x(t),t) = 0\). Differentiating the standard homotopy with respect to \(t\), we obtain \[0 = H'(x,t) = (1-t)g'(x)x' - g(x) + tf'(x)x' + f(x),\] so that we can describe \(x\) by \[x' = \frac{g(x) - f(x)}{(1-t)g'(x) + tf'(x)}.\]
As initial conditions, \(t=0\) and \(x\) any root of \(g(x)\).
For a general homotopy, including one in higher dimension, the differential equation takes the form \[x' = - \left(\frac{\partial H}{\partial x}\right)^{-1} \frac{\partial H}{\partial t},\] with the same initial conditions.
Example 12.3 Continuing Example 12.2, the associated differential equation is
\[x' = -\frac{x^2 + x - 12}{-3x^2 + (2t+4)x + (t+3)}\]
with any of the initial values \(x=-1,0,3\) at \(t=0\).
The most immediate concern for this differential equation would be places where the denominator is zero. If these poles cross our homotopy paths, we might be in trouble, so let’s take a look. This plot adds those poles as lines in red.
In this case, our root paths avoid any poles, so the differential equation should be fairly well-behaved.
12.3 Method
To find the root, we trace the path of roots by numerically solving the differential equation above. One of the great features of the method is that any numerical method for solving differential equations can be used. The version presented here is very simple, iteratively estimating with the BLA (Euler’s method) but it would not be difficult to substitute one of the better Unit 4 algorithms.
The BLA estimate is \[ x(t+\Delta t) \approx x(t) + x'(t)\Delta t, \]
which will be very accurate for a small \(\Delta t\). Starting from \(t_0=0\), we iteratively produce sequences \(t_n\) and \(x_n\). A natural choice is to run \(n\) steps with equally spaced steps, \(\Delta t = 1/N\). While \(x'\) also depends on \(x\), we can use our estimate \(x_n \approx x(t_n)\) in order to calculate \(x_{n+1}\)
Recall that this is exactly how your calculus classes arrived at the evenly-spaced left-handed Riemann sums – so we see that in the limit, it is exact, because \[\int_0^t x'(u) du = x(t) - x(0) = x(t)\] since \(x(0) = 0\) by assumption. So it’s not a new process for solving ODEs, the only twist is the added dependence on \(x\) at each step. When the step size is small, it produces passable estimates for most functions, without too much mathematical complexity.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
from shiny.express import input, render, ui
with ui.tags.head():
ui.tags.link(
rel="stylesheet",
href="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.css"
),
ui.tags.script(src="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.js")
ui.tags.style(
"""
body { background: #C8E0FF !important; }
.irs-bar { background: #0054A9 !important; }
.irs-handle { background: #0054A9 !important; }
.irs-single { background: #0054A9 !important; color: #FFFFFF !important; }
.irs-min { background:#99CBFF !important; color: #000000 !important; }
.irs-max { background:#99CBFF !important; color: #000000 !important; }
.control-label { display: block; text-align: center; }
.shiny-input-container { margin:auto; }
.shiny-input-select { background-color: #99CBFF; border: 1px solid #0054A9; }
.shiny-input-number { background-color: #99CBFF; border: 1px solid #0054A9; }
.shiny-input-number:focus { background-color: #99CBFF; border: 1px solid #0054A9; }
}
"""
)
def f(x):
return (x-2)*(x-3)*(x+2)
def g(x):
return x*(x+1)*(x-3)
def H(x,t):
return (1-t)*g(x) + t*f(x)
def dx(x,t):
return -((x**2) + x - 12)/(-3*x**2 + (2*t+4)*x + t + 3)
with ui.layout_columns(col_widths=(-1,4,-2,4,-1)):
ui.input_select("root_select", "initial root:", ["-1","0","3"])
ui.input_numeric("n_steps", "number of steps:", 10, min=1, max=100)
with ui.layout_columns(col_widths=(-2,4,4,-2)):
ui.p("estimated root:")
@render.ui
def est_root():
N = input.n_steps()
xn = int(input.root_select())
tn = 0
dt = 1/N
for i in range(0,N):
xn += dx(xn,tn)*dt
tn += dt
return ui.TagList(ui.tags.div(id="nm1"),ui.tags.script(f"katex.render(\"{xn:.4f}\", nm1);</script>"))In this example, with fairly nice functions, it still seems like quite a large number of steps are required to get decent accuracy; at \(100\) steps from \(x=-1\), only two digits after the decimal point are correct. The main cause is the accumulation of error from using the estimated \(x(t_n)\) when calculating \(x'(t_n)\). These snowball into a fairly large error across the whole interval. One way to start addressing this would be a more effective numerical method for solving the differential equation…
12.4 Path Correction
… but before jumping to new machinery, there’s an interesting and extremely efficient kind improvement to the method which leverages tools we’ve already developed.
Above, we solved the differential equation by a method that would work in great generality. In essence, we “forgot” that the differential equation came from \(H(x,t) = 0\), which gives it extra structure: we know that the actual solution comes from a root-finding problem. That is to say that at step \(n\) we have calculated some \(x_n\) so that \((x_n,t_n)\) is near the actual path, meaning \(H(x_n,t_n)\approx 0\). The true point on the path is the root of \(H(x,t_n)\). If we could adjust \(x_n\) to be closer to that root, then we’d not only make this particular point better, but also have better data for the next step of Euler’s method, reducing the amount of accumulating error from Euler’s method that we pointed out above (or any other method for solving the differential equation). This is path correction.
We already have a great method for refining a good guess to a better guess: Newton’s method. The path correction process is as follows:
- Carry out a step of Euler’s method to get a new \((x_n,t_n)\).
- Use Newton’s method to estimate a root of \(H(x,t_n)\) from the guess \(x_n\).
One could reasonably object that if we decide to spend computing power finding a root of \(H(x,t_n)\), we may as well just solve \(H(x,1)\) directly, perhaps by Newton’s method. However, we might not have a good guess for a root of \(H(x,1)\). In path correction, we know that Euler’s method should give us an \(x_n\) which is fairly close to the root of \(H(x,t_n)\), so it is extremly likely that that Newton’s method converges from this guess and does so rapidly. Just one or two steps of Newton’s method will increase accuracy significantly, and we’ll never have to come up with a guess ourself!
One step of Newton’s method requires only a little arithmetic, and the derivative of \(H(x,t_n)\) already has to have been calculated in order to write down the homotopy method’s differential equation; just evaluate \(\partial H/\partial x\) at \(t=t_n\).
12.5 Further Pictures
These two plots use our current example, plotting the plain (non-corrected) path using \(5\) and \(10\) steps of Euler’s method. The accuracy improves with 10 steps, but it’s not super accurate. In these examples, the actual root is \(2\).
Estimated root, 5 uncorrected steps:
Estimated root, 10 uncorrected steps:
Notice that in both uncorrected examples, the distance from the actual path increases at each step. This is typical.
With path correction, the results improve significantly. One path correction step roughly doubles the amount of calculation used, so 10 EM steps uncorrected is comparable to 5 EM steps corrected with 1 step of NM per EM.
Path-corrected root estimate with 5 steps of EM and 1 step of NM perm EM:
Path-corrected root estimate with 5 steps of EM and 2 step of NM per EM:
The path-corrected versions are so accurate that it’s not very easy to see both NM steps. Here’s a very coarse example with just three EM steps to make this more visible.
Path-corrected root estimate with 3 steps of EM and 2 step of NM per EM:
Notice that this is way more accurate than 10 steps of the uncorrected method!
The “actual” path in all of these examples was calculated by using 20 steps of Euler’s method with 2 path correction steps each :)
Exercises
Exercise 12.1 In our first example of numerically solving the homotopy continuation ODE, you may have noticed that the result for \(x=3\) was exact. Why?
Exercise 12.2 Let \(f(x)\) be some differentiable function and \(x_0\) a number. Define an infinite homotopy \[H(x,t) = f(x) - e^{-t}f(x_0)\] which has the property that \[\lim_{t\to\infty} H(x,t) = f(x)\] while \(H(x,0) = f(x)-f(x_0)\) has a root at \(x=x_0\). Therefore, we expect a path from \((x_0,0)\) to a root of \(f(x)\).
- Determine the differential equation for the homotopy method.
- Using Euler’s method with step size 1, write down a recursive formula which estimates the solution.
- What familiar numerical method is this?
Programming Problems
Exercise 12.3 Implement Euler’s method (evenly spaced) for implicitly defined functions, like the homotopy differential equation. Given \(\frac{dx}{dt} = f(x,t)\) and an initial point \((x_0,t_0)\), the Euler’s method estimate is \[x(t+\Delta t) = x(t) + f(x(t),t)\Delta t,\] which can be written as \[x_{n+1} = x_n + f(x_n,t_n) \Delta t\] where \(t_n = n\Delta t + t_0\).
Input:
- A function \(\frac {dx}{dt} = f(x,t)\).
- An initial point \((x_0,t_0)\) and end time \(t=b\).
- The number \(N\) of steps.
Output:
- An estimate for \(\int_{t_0}^b \frac{dx}{dt}(x,t)dt\).